import math

def xPrime(n):
    """
    Yield all prime up to n (n included)
    Basic sieve algorithm
    http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
    """
    # Check parameter
    if n < 2:
        return

    # Handle even prime number
    yield 2

    # Create a odd number prime property
    #   value     1,  3,  5,  7,  9, 11, 13, 15, ...
    #   sieve = [ F,  T,  T,  T,  F,  T,  T,  F, ...]
    #
    # Offset manipulation
    #   p       = 2*x + 1            == sieve_offset(x)
    #   p^2     = 4*x^2 + 4*x + 1
    #           = 2(2*x*(x + 1)) + 1 == sieve_offset(2*x*(x + 1))
    #   p + 2*a = 2*x + 1 + 2*a
    #           = 2*(x + a) + 1      == sieve_offset(x + a)
    size  = (n + 1) // 2
    sieve = [True] * size

    # Compute sieve size limit
    root = (int(math.sqrt(n)) + 1) // 2

    # Sieve and return odd prime number
    for i in xrange(1, root):
        if sieve[i]:
            # sieve[i] = 2.i + 1 = p is prime
            yield 2*i + 1

            # remove odd multiple (p^2 + 2.k.p for k=0...)
            # p^2 = 4.i.i + 4.i + 1 = 2.2.i.(i + 1) + 1
            # p^2 = sieve[2.i.(i+1)]
            # sieve[x] + 2.p = s[x + p] = s[x + 2.i + 1]
            for j in xrange(2*i*(i + 1), size, 2*i + 1):
                sieve[j] = False

    for i in xrange(root, size):
        if sieve[i]:
            yield 2*i + 1


def Prime(n):
    """
    List all prime up to n (n included)
    """
    return [p for p in xPrime(n)]


def _TryComposite(a, d, s, n):
    """
    Check composite for Miller-Rabin primality test
    Retrun True if n is composite
    http://en.wikipedia.org/wiki/Miller-Rabin_primality_test
    """
    x = pow(a, d, n)
    if x == 1 or x == n - 1:
        return False
    for i in xrange(s - 1):
        x = pow(x, 2, n)
        if x == 1    : return True
        if x == n - 1: return False
    return True

def IsPrime(n, precision=10):
    """
    Miller-Rabin primality test
    Return True if n is prime otherwise False
    http://en.wikipedia.org/wiki/Miller-Rabin_primality_test
    """
    # Simple case
    if n < 2:
        return False
    if n in [2, 3, 5, 7, 11, 13, 17, 19]:
        return True
    if any(n % p == 0 for p in [2, 3, 5, 7, 11, 13, 17, 19]):
        return False

    # Write n-1 as d * 2^s wth d odd
    d, s = n - 1, 0
    while d % 2 == 0:
        d, s = d >> 1, s + 1

    # try composite
    if   n <         1373653:   al = [2, 3]
    elif n <         9080191:   al = [31, 73]
    elif n <      4759123141:   al = [2, 7, 61]
    elif n <   1122004669633:   al = [2, 12, 23, 1662803]
    elif n <   2152302898747:   al = [2, 3, 5, 7, 11]
    elif n <   3474749660383:   al = [2, 3, 5, 7, 11, 13]
    elif n < 341550071728321:   al = [2, 3, 5, 7, 11, 13, 17]
    else:                       al = Prime(100)[:precision]
    return not any(_TryComposite(a, d, s, n) for a in al)

